library(list)
library(ggplot2)

###load March dataset
load("MarCourier.RData")

#################Variable information
####All variables with a political figure's name indicate support for the figure (1=support, 0=not support)

####Contemporary is the number of contemporary figures supported on a list, History the number of 
#historical figures supported on a list, International the number of international figures supported in a 
#list (the Castro experiment), and Placebo the number of items relevant to an individual (the Placebo 
#experiment)

###TContemporary, THistory, TPlacebo and TInternational are vectors representing whether or not an individual 
#was in a control or treatment group for the respective experiments (1=Treatment, 0=Control)

###Age is a vector of the respondents' age

###################Analyses
##########################################
########################Analysis of Placebo, Table 2
##########################################
t.test(mar$Placebo[mar$TPlacebo==0],mar$Placebo[mar$TPlacebo==1])
##Compare to number of individuals who reported being over 55 in direct qs
mean(mar$Age>55)

###Blair Imai test of design effect; transform into vectors for condition and results, remove nas, fn 25

pna<-mar$Placebo[is.na(mar$Placebo)==F]
tna<-mar$TPlacebo[is.na(mar$Placebo)==F]

ict.test(pna,tna,3)

####frequency table, Table A1
table(mar$Placebo,mar$TPlacebo)


################################################
###################Analysis of Castro experiment, table 2
################################################
t.test(mar$International[mar$TInternational==0],mar$International[mar$TInternational==1])
##compare to direct support for Castro
mean(na.omit(mar$Castro))

###Blair and Imai test for design effects, same technique as in placebo, fn 25
ina<-mar$International[is.na(mar$International)==F]
tna<-mar$TInternational[is.na(mar$International)==F]
ict.test(ina,tna,3)

####frequency table, Table A1
table(mar$International,mar$TInternational)

###Analyze number of figures supported in control group in direct questions and list
####vector os support for all figures
itf<-mar$Mandela+mar$Lukashenko+mar$Merkel

####Analyze differences between control population in terms of direct and list support, fn 26
mean(na.omit(itf[which(is.na(mar$International)==F & mar$TInternational==0)]))
mean(na.omit(mar$International[which(is.na(itf)==F& mar$TInternational==0)]))

#############################################
####################Putin analyses
#############################################

##support for Putin in direct question, Table 1
mean(na.omit(mar$Putin))

##################################
###########Historical analysis, Table 1
#################################
t.test(mar$History[mar$THistory == 0],mar$History[mar$THistory == 1])


####Test for design effect, create vectors for experiment results and condition, remove nas, run Blair and Imai, fn 25 
hna<-mar$History[is.na(mar$History)==F]
tna<-mar$THistory[is.na(mar$History)==F]
ict.test(hna,tna,3)


####frequency table, Table A1
table(mar$History,mar$THistory)

###vector of support for all figures
mar$hf<-mar$Stalin + mar$Brezhnev + mar$Yeltsin


################Floor effects
########Support for Putin vs. all others, not reported
mean(na.omit(mar$Putin == 1 & mar$Stalin == 0 & mar$Brezhnev == 0 & mar$Yeltsin==0))

###percent at 1, treatment, historica, fn 22
mean(na.omit(mar$History[mar$THistory==1])==1)

###modify historical experiment so that individuals at one in treatment receive 0; fn 22
H2<-mar$History
H2[which(mar$History==1 & mar$THistory==1)]<-0
t.test(H2[mar$THistory==0],H2[mar$THistory==1])

### difference between individuals in control in number supported, nas removed; experiment first; not reported, akin to fn 26
mean(na.omit(mar$History[which(is.na(mar$hf)==F & mar$THistory==0)]))
mean(na.omit(mar$hf[which(is.na(mar$History)==F & mar$THistory==0)]))

###deflation regression, Table A2
summary(lm(History ~ THistory*hf, data=mar))

####Analyze relationship betweeen number of control figures supported and experimental results, Figure 2
Condition<-as.factor(mar$THistory)
levels(Condition)<-c("Control", "Treatment")
qplot(mar$hf,mar$History, xlab="Number of control-group figures supported in direct questions", 
      ylab="Number of figures supported in list", alpha=I(.0001)) +
  geom_smooth(aes(linetype=Condition), color="black", size=1) + theme_bw() + theme(legend.position=c(.2,.8)) +
  scale_linetype_manual(values=c("solid","dotted"))

###Get predictions of SDB, following Blair and Imai
set.seed(123)
summary(dim.r <- ictreg(History ~ 1, treat = "THistory", data=mar,J=3, method = "lm"))
pdim.r<-predict.ictreg(dim.r)

df<- glm(Putin ~ 1, data=mar,  family=binomial)

predict(dim.r, direct.glm=df)

#################################
###############Contemporary analysis, table 1
###########################
t.test(mar$Contemporary[mar$TContemporary == 0],mar$Contemporary[mar$TContemporary == 1])



####Test for design effect, create vectors for experiment results and condition, remove nas, run Blair and Imai analysis, fn 25

cna<-mar$Contemporary[is.na(mar$Contemporary)==F]
tna<-mar$TContemporary[is.na(mar$Contemporary)==F]
ict.test(cna,tna,3)


####frequency table, Table A1
table(mar$Contemporary,mar$TContemporary)

###vector of support for figures
mar$cf<-mar$Zyuganov + mar$Zhirinovskii + mar$Mironov

#######Floor analysis
######## Putin vs. all others, not reported, akin to fn 15
mean(na.omit(mar$Putin == 1 & mar$Zhirinovskii == 0 & mar$Zyuganov == 0 & mar$Mironov==0))

###percent at 1, treatment, contemporary, fn 22
mean(na.omit(mar$Contemporary[mar$TContemporary==1])==1)

###modify contemporary experiment so that individuals at one in treatment receive 0, fn 22
C2<-mar$Contemporary
C2[which(mar$Contemporary==1 & mar$TContemporary==1)]<-0
t.test(C2[mar$TContemporary==0],C2[mar$TContemporary==1])

 
###determine difference between control number of figures supported in list and direct questions, fn 26
mean(na.omit(mar$Contemporary[which(is.na(mar$hf)==F & mar$TContemporary==0)]))
mean(na.omit(mar$cf[which(is.na(mar$Contemporary)==F & mar$TContemporary==0)]))

###deflation regressions, Table A2
summary(lm(Contemporary ~ TContemporary*cf, data=mar))

####Analyze relationship betweeen number of control figures supported and experimental results, Figure 2
Condition<-as.factor(mar$TContemporary)
levels(Condition)<-c("Control","Treatment")

qplot(mar$cf,mar$Contemporary, xlab="Number of control-group figures supported in direct questions", 
      ylab="Number of figures supported in list", alpha=I(.0001)) +
  geom_smooth(aes(linetype=Condition), color="black", size=1) + theme_bw() + theme(legend.position=c(.2,.8)) +
  scale_linetype_manual(values=c("solid","dotted"))

###Get predictions of SDB with Blair/Imai method, Table 1

set.seed(123)
summary(dim.r <- ictreg(Contemporary ~ 1, treat = "TContemporary", data=mar,J=3, method = "lm"))
pdim.r<-predict.ictreg(dim.r)

df<- glm(Putin ~ 1, data=mar,  family=binomial)

predict(dim.r, direct.glm=df)

##################################
##############x2 list analysis, Table 1
#################################

A<-cbind(mar$History[mar$TContemporary==1],mar$Contemporary[mar$TContemporary==1]) ##Create matrix for treatment
B<-cbind(mar$History[mar$TContemporary==0],mar$Contemporary[mar$TContemporary==0]) ##create matix for control

##lw remove nas
A<-A[!is.na(A[,1]), ]
A<-A[!is.na(A[,2]), ]
B<-B[!is.na(B[,1]), ]
B<-B[!is.na(B[,2]), ]

###CIs from Glynn method
(mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2

(mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2+qnorm(.975)*sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))
(mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2-qnorm(.975)*sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

